Non-equilibrium cluster-perturbation theory 
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The cluster perturbation theory (CPT) is one of the simplest but systematic quantum cluster 
approaches to lattice models of strongly correlated electrons with local interactions. By treating the 
inter-cluster potential, in addition to the interactions, as a perturbation, it is shown that the CPT 
can be reformulated as an all-order re-summation of diagrams within standard weak-coupling per- 
turbation theory where vertex corrections are neglected. This reformulation is shown to allow for a 
straightforward generalization of the CPT to the general non-equilibrium case using contour-ordered 
Green's functions. Solving the resulting generalized CPT equation on the discretized Keldysh- 
Matsubara time contour, the transient dynamics of an essentially arbitrary initial pure or mixed 
state can be traced. In this way, the time-dependent expectation values of one-particle observables 
can be obtained within an approximation that neglects spatial correlations beyond the extension of 
the reference cluster. The necessary computational effort is very moderate. A detailed discussion 
and simple test calculations are presented to demonstrate the strengths and the shortcomings of the 
proposed approach. The non-equilibrium CPT is systematic and is controlled in principle by the 
inverse cluster size. It interpolates between the non-interacting and the atomic or decoupled-cluster 
limit which are recovered exactly and is found to predict the correct dynamics at very short times 
in a general non-trivial case. The effects of initial-state correlations on the subsequent dynamics 
and the necessity to extend the Keldysh contour by the imaginary Matsubara branch are analyzed 
carefully and demonstrated numerically. It is furthermore shown that the approach can describe 
the dissipation of spin and charge to an uncorrelated bath with an essentially arbitrary number of 
degrees of freedom. 

PACS numbers: 71.10.Fd, 71.27. +a, 67.85.Lm 



I. INTRODUCTION 



A theoretical understanding of transient processes in 
systems of strongly correlated electrons far away from 
thermal equilibrium and the development of accord- 
ing methods is one of the most challenging tasks in 
condensed-matter physics. There is a pure theoretical 
motivation, on the one hand, since the study of non- 
equilibrium states opens up a new perspective on clas- 
sical many-body effects, such as collective magnetic or- 
der, high-temperature superconductivity, Kondo screen- 
ing of local moments or Mott metal-insulator transi- 
tions, for example. On the other hand, there is an 
urgent need to describe and understand the results of 
recent exciting experimental studies in different fields: 
This includes nanostructure physics as, for example, 
the application of scanning-tunnelling microscope tech- 
niques to measure the spin relaxation time of itinerant 
and correlated electrons in nanostructures with atomic 
resolution^ or relaxation and switching times in first 
atom-by-atom realizations of all-spin based spintron- 
ics devices^ Furthermore, an improved theoretical un- 
derstanding of fast demagnetization processes probed 
by femtosecond optical excitations^ and of the non- 
equilibrium electronic structure of strongly correlated 
transition-metal oxides which is accessible to femtosec- 
ond pump-probe spectroscopies^^ Another fascinating 
field is the controlled preparation and monitoring of 
the non-equilibrium dynamics of highly excited fermionic 
states realized in correlated systems of ultracold atoms 
in optical lattices^ In all these examples, the most in- 



teresting questions refer to the effect of strong nonlocal 
electronic correlations on the dynamics of itinerant elec- 
trons on a lattice or in a well-defined nanostructure. 

For the strong-correlation regime of extended sys- 
tems, non-perturbativc numerical methods are re- 
quired. Besides exact-diagonalization techniques^ which 
are limited to systems with small Hilbcrt-space di- 
mensions, numerical renormalization-group^ or density- 
matrix renormalization-group techniques 9 - can be used to 
study impurity or one-dimensional lattice systems with 
high numerical accuracy. The continuous-time quan- 
tum Monte-Carlo approach can straightforwardly be ex- 
tended to the non-equilibrium caseJ^ It also belongs to 
the class of numerically exact methods but is limited, 
due to the sign or phase problem, to short-time dynam- 
ics. Among the non-perturbative but approximate tech- 
niques, Green's function-based embedding methods are 
attractive. Relying on the pioneering work of Kubo^ 
Schwinger^ Kadanoff, BaymI2 and Keldysh^ (see also 
Refs. Il5l4l7l) all-order diagrammatic re-summations can 
be used to define non-equilibrium generalizations of 
dynamical mean-field theory ; 18 ' 19 self-energy functional 
theory^ or of the dual-fermion approach^ All of the 
above-mentioned impurity or cluster-embedding methods 
are highly expensive numerically. 

The purpose of the present paper is to propose and 
to discuss a method which is obtained by a generaliza- 
tion of the cluster-perturbation theory (CPT)»22~— This 
non-equilibrium CPT is a conceptually simple method 
which can be applied to lattice models of correlated elec- 
trons with local interactions and basically arbitrary ini- 
tial states and arbitrary Hamiltonian dynamics. The re- 
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quired computational resources are very moderate. It 
is based on a partitioning of the lattice model of in- 
terest into smaller parts ("clusters") that are amenable 
to an exact solution, preferably by means of exact di- 
agonalization, and treats the initially disregarded inter- 
cluster terms subsequently in an approximative way such 
that the method becomes systematic and controlled by 
the inverse cluster size. The non-equilibrium CPT ac- 
counts for temporal correlations and includes non-local 
but short-range spatial correlations up to the scale of the 
cluster size in the spirit of cluster mean-field methods^ 
It is thereby closely related to the (cellular) dynamical 
mean-field approach, and can be seen as the starting 
point for more elaborate but also more expensive self- 
cncrgy-functional or dual-fermion techniques. The pro- 
posed non-equilibrium CPT is the simplest systematic 
approach to non-equilibrium dynamics which includes 
non-local correlations. 

Our formal idea is to first re-construct the usual equi- 
librium CPT by means of the standard weak-coupling 
perturbation expansion but treating besides the bilinear 
inter-cluster hopping the quartic interaction terms as a 
perturbation as well. The CPT Green's function is then 
obtained by formally summing all diagrams to infinite 
order but neglecting certain vertex corrections. In a sec- 
ond step, this idea can straightforwardly be transferred 
to the non-equilibrium situation by replacing the ther- 
mal Green's function with the contour-ordered Green's 
function. The central CPT equation thereby becomes a 
matrix equation in orbital and (discretized) time indices 
which can easily be solved numerically. 

The paper is organized as follows: The basic theory of 
non-equilibrium Green's functions is reviewed in the next 
section |TT] with notations following Ref. [l?! Section IIIII 
develops the non-equilibrium cluster-perturbation theory 
in detail. An extensive discussion of the new approach 
and of different numerical results is given in section IIVI 
The conclusions arc summarized in section IVl 



II. EXPANSION OF THE NON-EQUILIBRIUM 
GREEN'S FUNCTION 



Consider a system of electrons which at time to is in 
a normalized pure state |^). We assume that this state 
is the TV-particle ground state of some properly defined 
Hamiltonian 



B — -Bo + B\ , 



(1) 



where Bq is a one-particle operator and B\ an interaction 
term. Alternatively, the system could be at time to in a 
mixed state p where it is assumed that a Hamiltonian B 
can be found such that 

P trexp(-/3£) ' { 1 

where B = Bo + B\ = B — p,N and where /3 is the inverse 
temperature of the initial state. With p = and 
(3 — > oo this also comprises pure initial states. 
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FIG. 1: Three-branch contour 7 in the complex time plane. 
The upper and the lower real branches define the Keldysh con- 
tour, the imaginary branch is called the Matsubara branch. 



For t > to the system's time evolution shall be gov- 
erned by the explicitly time-dependent Hamiltonian 



H(t) = H +H 1 (t), 



(3) 



where [H(t),B]- ^ in general. For the calculations 
below, we will assume that the system is not driven by 
explicitly time-dependent external fields and that H (t) = 
H(0) ^ B. However, the formalism will be developed for 
the general case. 

Consider an arbitrary possibly time-dependent observ- 
able A(t). Its time dependence within the Heisenberg 
picture with respect to H(t) = H(t) — pN is determined 
by the equation of motion 

ij t Au(t) = [A H (t), H(t)}- + i^A H (t) (4) 

with the initial condition A-u(to) = A(to). The formal 
solution of the equation of motion is given by 



A n (t) = [Te^o^l A(t) \ Te -if: *™® 



(5) 



where T (T) is the chronological (anti-chronological) 
time-ordering operator. 

For a system in the initial state p the expectation value 
of the observable A(t) at time t is (A) t = tr(pA-u(t)). 
This can be written &S)H 



(A)t = 



tr (% cxp (-i J 7 dt JC(t)) A{t) 
tr (r y cxp (-i / n 



dtK{t) 



(6) 



Here, the time integration is carried out along the contour 
7 in the complex time plane. 7 extends from t = to to 
t = 00 along the real axis (upper branch) and back from 
t = 00 to t = to along the real axis (lower branch) and 
finally from t = to to t = to — i/3 along the imaginary 
axis (Matsubara branch), see Fig. [T] We also refer to 
the upper and the lower branch as the Keldysh contour. 
%f denotes the ordering operator along the contour and, 
after expanding the exponential, places an operator JC(t\) 
to the left of ICfa) if t\ is "later" than t% on the contour 7 
where to ~ i/3 is the "latest" time. Obviously, T n replaces 
T on the upper and T on the lower branch. Finally, 
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K,(t) = H{t) on the upper and the lower branch of 7 
while K(t) — B on the Matsubara branch. 

T-y also acts on A(t). The time argument of A{t) 
is the time at which the expectation value is evalu- 
ated and indicates the position on the time contour 
where for the integrals in the numerator the observ- 
able has to be placed. Note that, for the numerator, 
the results of integrating along the upper and the lower 
branches between t = t and t = 00 cancel each other. 
Hence, the integration along the Keldysh part of the 
contour can be limited to t < t. For the denominator, 
only the Matsubara branch of the contour contributes 
to the integral with the result trexp(— (3B). If H is 
time-independent and equal to B, the equilibrium result 
(A) t = tr{exp(-PH)A(t))/trexp(-PH) is recovered. 

We assume the Hamiltonian B, which characterizes the 
initial state, and the Hamiltonian H(t), which determines 
the system's dynamics, to be given in second-quantized 
form: 

B = E T S ] ^ + \ E u ZA4 c ^ ( 7 ) 

ex/3 a/3j8 

and 

# (*) = E T »p c i°0 + 2 E 0a/w 7 (t)44 G r < * ■ ( 8 ) 

a/3 a/3~f5 

Here a refers to a complete and orthonormal set of (time- 
independent) one-particle orbitals, i.e. the explicit time- 
dependence is due to the interaction parameters only. An 
external bilinear time-dependent field could be consid- 
ered in addition. In this case the interaction part would 
also contain terms bilinear in c a and . 

The time-dependent expectation value of any one- 
particle observable A(t) = J2 a aa l3(^) c a c P can ^ e 
tained from the contour-ordered Green's function 

iG aa ,{t,t') = (%c K , a {t)J Kia ,{t')) (9) 

as 

(A)t = -iJ2a<xt3(t)Gp a (t,t + + ) (10) 

a/3 

where + is a positive infinitesimal and (• • • ) = tr(p • • • ) 
denotes the expectation value in the initial state. Fur- 
thermore, the annihilator and the creator are given in 
the Heisenberg picture with respect to Kit), t,t' are ar- 
bitrary times on the contour, and 7^ is the time ordering 
of annihilators and creators on the contour 7 which yields 
an additional (Fermi) sign per transposition. 

The contour-ordered Green's function involves opera- 
tors given in the Heisenberg picture, i.e. with a time- 
dependence due to the interacting Hamiltonian TL, and 
an expectation value with a (mixed) state corresponding 
to the interacting Hamiltonian B. The main motivation 
for placing the contour-ordered Green's function in the 
focus of the theory, rather than, for example, expecta- 
tion values like (A) t , is that (i) the Green's function can 



be brought into a form that meets the requirements to 
apply Wick's theorem and that (ii) the application of 
Wick's theorem only generates contour-ordered Green's 
functions again. Thereby, a closed set of physically in- 
teresting quantities is obtained, and a consistent pertur- 
bation theory can be set up. 

Following Ref. |T3, the contour-ordered Green's func- 
tion can be cast into the form: 

. n u _ (T, e- i Afeo- 1 ^ CCo , a (t)4 0|a ,(tQ)W 

i(j aa i (t,t ) — . r ,~* • 

(T e~ l ^i dt?c ' c o.iW)(o) 

(11) 

In this expression, the annihilators and creators, cje 0)O ,(t) 
and C£ a '(t') possess a "free" time dependence only, i.e. 
they are given in the interaction picture where the time 
dependence is due to /Co only. The same applies to the in- 
teraction term /C/c 0i i(i) appearing under the contour inte- 
gral - its time dependence is "free" and given by /Co only. 
Finally, also the expectation value (• • • )(°) = tr(po ■ ■ • ) is 
a "free" one and is defined with free density operator 
po = cxp(—f3Bo)/Z only. Hence, we can apply Wick's 
theorem and therewith standard techniques of perturba- 
tion theory 

Expanding the Green's function in powers of the inter- 
action parts of B and H, the n-th order coefficient turns 
out to be given in terms of 2n + 1 "free" contour-ordered 
Green's functions: 

i(^M = (TyC^At)^,^)) l0) - (I 2 ) 

This can be computed exactly for the case considered 
here, i.e. for U (t) = H = const, but [B ,H ]- + 0. We 
find: 

aqW (t +>\ - ( r -i(T K -p,)t 1 i(T K -u,)t' ) 

^aa'K Z ^)-\ e 1 + e -P(T B -n) e J 

\ / CLOi' 

(13) 

if t later than t' on 7 and 

7(7 (0) (t t') - - ( r -i(T K -p.)t 1 i(T K -p,)t'\ 

\ e^B-M) + i e ) aa , 

(14) 

if t' later than t on 7. On the Keldysh contour t is real 
and Tk = T with the elements T a p while on the Matsub- 
ara branch t = —it with < r < /3 and Tk = Tb with 
elements see Eq. © and Eq. ©. 

III. CLUSTER-PERTURBATION THEORY 

There are several ways to define the cluster- 
perturbation theory (CPT) for the equilibrium case. 
The first approach, based on the so-called Hubbard-I 
approximation^ focuses on the electron self-energy of 
the Hubbard model2£~— for a D dimensional lattice. The 
Hubbard-I approximation can be constructed by starting 
from the atomic limit of the Hubbard model and taking 
the self-energy from that limit as an approximation for 
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the infinite lattice model. In the original work^ addi- 
tional requirements on the average occupation numbers 
are imposed which must be solved sclf-consistcntly. The 
Hubbard-I concept was generalized later— by starting 
from a finite Hubbard cluster instead of a single Hub- 
bard atom. Approximating the lattice self-energy by the 
self-energy of a cluster consisting of a finite number of 
L c sites, defines a numerical technique which (i) directly 
works in the thermodynamical limit, (ii) can be improved 
systematically by increasing the cluster size L c and (iii) 
provides, via Dyson's equation, a single-electron Green's 
function which respects certain general requirements of 
Lehmann representability and causality. On the other 
hand, this construction of the CPT appears to be rather 
ad hoc. 

The second approach is based on strong-coupling per- 
turbation theory for the Hubbard model and is more sys- 
tematic. For Hubbard-type models, an expansion in pow- 
ers of the hopping t around the atomic limit can be or- 
ganized in a systematic diagrammatic series <2£r— At the 
lowest order this leads to the Hubbard-I approximation. 
The CPT is obtained from a cluster generalization of the 
strong-coupling expansion. The extension consists in a 
partitioning of the lattice into small clusters that can be 
treated exactly, and a subsequent expansion in powers of 
the inter-cluster hopping. The lowest order constitutes 
the CPTi 22 ' 23 In principle, the expansion can be carried 
out to arbitrary order in the inter-cluster hopping using 
the diagrammatic method of Refs. I30f33l or the cluster 
dual-fermion method^ However, going beyond the low- 
est order is quite demanding numerically and leads to 
causality problems at large t and low temperatures due 
to the degeneracy of the ground state. Since the low- 
est order of the strong-coupling expansion is causal and 
still represents a systematic approach with respect to the 
cluster size L c , it has gained some attraction in the past. 
The CPT is a conceptually simple method which never- 
theless includes short-range correlations on the scale of 
the cluster size and which requires moderate computa- 
tional resources only. 

An alternative approach to construct the CPT is pre- 
sented here. It is based on the usual weak-coupling 
perturbation expansion. Besides the quartic interaction 
terms in B and however, we additionally treat the 

bilinear inter-cluster hopping perturbation as well. 
The CPT Green's function is then obtained by formally 
summing all diagrams to infinite order but neglecting 
vertex corrections. This idea can straightforwardly be 
transferred to the non-equilibrium situation by replacing 
the thermal Green's function with the contour-ordered 
Green's function. 

Starting point for the construction of the CPT is a par- 
titioning of the original D dimensional lattice consisting 
of L sites into clusters of finite size and open bound- 
aries. The clusters shall consist of L c sites each. Fig. 
Fig. [U gives an example for the D = 2 square lattice and 
L c = 4. For simplicity, we assume all clusters to be iden- 
tical and to form a supcrlatticc labeled by a superlattice 
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FIG. 2: Partitioning of a D = 2 square lattice into clusters 
with L c = 4 sites each. Ti : ij denotes the intra-cluster hopping 
between sites i and j within the same cluster /. Vjj.y is the 
inter-cluster hopping between sites i £ I and j € J. 

site index 7=1,..., L/L c . The sites within the cluster I 
are labeled by an index i = 1, L c . 

The Hamiltonians of the initial thermal and of the 
transient final state, i.e. B and H(t), are decomposed 
accordingly, 

B = B' + B inter , H(t)=H'(t)+H intcI . (15) 

B' and H'(t) correspond to the reference system of de- 
coupled clusters. We have 

L/L c L/L c 

B' = £ B'j , H'{t) = £ H'M , (16) 
i=i i=i 

where B\ and H'j(t) describe the thermal initial state 
and the dynamics of the isolated cluster I. The CPT is 
mainly designed for applications to Hubbard-type models 
with local interactions. Besides the intra-cluster hopping, 
we therefore assume the interaction terms B\ and H\ (t) 
to be fully included in the reference system. Hence: 

B 'i=J2 X • (17) 

and 

l c 

H 'i® = EE £ /,iM^ c L 4 c««r, + Hij(t) , (18) 

i,j = l <T;<Tj 

where i and j run over the sites within the cluster I, 
where <Tj labels the residual orbital and spin degrees of 
freedom at a site i, and where B\j and H\j(t) denote 
the respective interaction part within cluster I. On the 
other hand, the inter-cluster parts include bilinear hop- 
ping terms only: 

-Sinter = Mntcr,IJ 7 -fainter = / ] Writer, I J (19) 

I..J I,J 

where 

Bin t er,IJ= £ V uk^A^ ' ^ 

ieijeJ en <jj 
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and 

H inter, I.I = ^2 ^2 \ - <t '\, r <> . (21) 

A triple of indices (I,i,ai) = a labels a certain or- 
bital of the one-particle basis. With respect to this ba- 
sis, the intra-cluster and the inter-cluster hopping pa- 
rameters form matrices £b, Vb and e, V, respectively. 
We have T B = e B + V B and T = e + V, see Eq. © 
and Eq. ([8]). In case that the superlattice of clusters 
is invariant under translations, Fourier transformation 
block-diagonalizes Vb and V simultaneously. Exploiting 
the fact that the intra-cluster hopping is already diag- 
onal in and independent of the superlattice index /, we 
get matrices of the form: eb, Vb(&) and e, V(fc), respec- 
tively, where k has the physical meaning of a wave vector 
and where sb and e are wave- vector independent. In all 
other cases, diagonalization must be done numerically, if 
desired. Note that Tb and T are different for a general 
initial state and cannot be diagonalizcd simultaneously. 

To set up the perturbation theory based on Wick's the- 
orem, the quartic terms Bi and Hi (i) have to be treated 
as a perturbation. As concerns the bilinear terms £>o and 
Hq, however, we are free to treat them as "free" or as 
a "perturbation" . Any choice is consistent with Wick's 
theorem. A non-equilibrium generalization of the CPT 
is obtained when treating the inter-cluster couplings Vb 
and V as perturbations while Sb and e are considered 
to be free. 

Perturbation theory then provides us with Dyson's 
equation for the fully interacting contour-ordered Green's 
function: 

G = G' + G' .H UK , Vk [G' ]-G. (22) 

Here, all quantities are matrices with respect to time vari- 
ables and orbital indices, such that the Green's function 
G has the elements G aa > (t, t'), for example, and Eq. (|22|) 
is short for: 

G aa ,{t,t') = G^J(t,t')+ I I dt " dt "' 

Ol"Ol"' * 1 

(23) 

The free Green's function G' in Eq. ([22]) is the U B = 
U(t) = 0, Vb = V = Green's function, i.e. 
the interaction-free intra-cluster contour-ordered Green's 
function or the interaction-free Green's function of the 
reference system. Explicitly, we have: 

j A i + e -^B-M) e ) aa , 

(24) 

if t later than t' on 7 and 

iG {0) '(t,t') = -( e-^-^ 1 —, — e <(*K-M)A 

act y e/K.B-fO+1 J m , 

(25) 



if t' later than t on 7. Here, £k = £b if t = —it is 
on the Matsubara branch and £k = £ for real t on the 
Kcldysh contour. The self-energy E^ Ki v K [G' ] inEq. (|22|) 
is obtained by summing over all irreducible self-energy 
insertions, formed by free propagators G' and vertices 
U K (t) and V K where U K (t) = U B or U K {t) = U{t), 
and likewise for Vk, depending on the position of the 
respective vertex on the time contour. 

The exact self-energy can formally be obtained in a 
two-step renormalization procedure, see Fig. [3Jt. First, 
we consider the renormalization of the free propagators 
due to Vk, i.e. due to electron scattering at the non- 
local but instantaneous (local in time) inter-cluster po- 
tential. The corresponding self-energy is simply given 
by Sv K [G' ] = Vk <8> 1 with the <5-function on the con- 
tour l tjt < = <5 7 (t, t'), and the renormalized propagator Go 
is obtained as the solution of the corresponding Dyson 
equation: 

Go = G + G • V K ® 1 • G* . (26) 

This yields the Green's function for 17b = U(t) = 0. Sub- 
sequent J7k renormalization is formally achieved by intro- 
ducing the corresponding self-energy S[/ K [^0] which is a 
(highly complicated) functional of the VK-renormalized 
propagator. This yields the full propagator as the solu- 
tion of 

G = G + Go-S [/K [Go]-G. (27) 

Since all diagrams are summed up, the procedure is ex- 
act. Comparison with Dyson's equation Eq. (|22[) shows 
that 

Sf/K,v K [Go] = S [/k [(G " 1 -Vk®1)- 1 ] + V k ®1. (28) 

Since the self-energy S[/ K [Go] is essentially unknown, 
this does not provide, of course, a pragmatic way to com- 
pute the full propagator. 

Let us now consider the Uk renormalization first, see 
Fig. \3]p. This leads to the following Dyson equation: 

G' = Go + G -5WGo]-G'. (29) 

Its solution is the interacting Green's function G' of the 
reference system of decoupled clusters. While still the 
functional form of S[/ K [G ] is highly non-linear and un- 
known, the propagator G' may be calculated directly by 
exact diagonalization, provided that the cluster size L c 
is moderate. Note that here it is essential to assume the 
J7k vertex to be local and not to couple different clusters. 
The subsequent Vk renormalization of the already Uk~ 
renormalized propagators is expressed with the Dyson 
equation 

Gcpt = G' + G' ■ V K ® 1 ■ G CPT . (30) 

Its solution defines the non-equilibrium CPT Green's 
function Gcpt- The reversed two-step renormalization 
is not exact since there are certain diagrams missing, see 
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FIG. 3: Re-summation of diagrams generated by electron- 
electron scattering U and by scattering at the inter-cluster 
potential V. (a) Exact procedure: Renormalization of the 
free (U = V = 0) propagator G' by potential scattering fol- 
lowed by renormalization of the U = propagator Go due to 
electron scattering [Eqs. (T2(j]) and [[27))]. (b) CPT: Renormal- 
ization of the free (U = V = 0) propagator G'o by electron 
scattering followed by renormalization of the V = propaga- 
tor G' due to potential scattering [Eqs. ()29|) and (|30|l ]. (c) 
Self-energy diagram, second order in U, second order in V, 
which is not taken into account within CPT. 



Fig. [3J;. From Eq. 
sclf-cncrgy 



and Eq. ([H]) we get the CPT 



Scpt^k.VkIGq] = X Uk [G' ] + V K ® 1 



(31) 



Comparing this expression with the exact self-energy Eq. 
(f2"5| shows that CPT neglects the influence of scattering 
at the inter-cluster potential on the renormalization of 
propagators due to the interaction, i.e. vertex correc- 
tions. Another way to paraphrase the approximation 
is to say that the CPT neglects electron-electron (17k) 
scattering across different clusters but takes into account 
intra-cluster electron-electron scattering and scattering 
of electrons dressed by Uk processes at the one-particle 
inter-cluster potential. 



IV. DISCUSSION AND RESULTS 

In the following we discuss the non-equilibrium CPT 
in detail and present numerical results to demonstrate, 
as a proof of principle, that the approach can be used in 
practice. 



A. Thermal equilibrium 

First, it has to be shown that the usual CPT is recov- 
ered for the case of thcrmodynamical equilibrium. We 
therefore assume that H(t) = H = B for a moment. 
Inspection of Eqs. (J23J) and ([2SJ and of Eqs. (Q2]) and 



P| immediately shows that G^, (t, t') and G^, (t, t') 
become temporally homogeneous, i.e. become functions 
of t — t' only. The interacting Green's function of the 
reference system, G' aa ,(t,t'), has to be computed exactly 
within non-equilibrium CPT and, therefore, is homoge- 
neous. Since Vk ® 1 is homogeneous by definition, the 
CPT equation (Eq. (|3"0J)) proves the CPT Green's func- 
tion Gcpt to be homogeneous, too. With Eq. (|T0|) this 
implies that the expectation value of any (not explicitly 
time-dependent) observable A is constant, (A) t = (A) to , 
and given by its thermal value for all t > U). 

Furthermore, as is shown below, there is an indepen- 
dent CPT equation on the Matsubara branch only: 



*2cpt — G 



G' ■ V B (8 1 • G 



1j +Lx • VB 



A iiCPT 



(32) 

Here the underlined symbols represent matrices in t, t' 
(besides orbital indices) where t, t' are restricted to the 
Matsubara branch only and where the integrations im- 
plicit in the notations are limited accordingly. Together 
with the homogeneity of the quantities, this allows to 
transform to a Matsubara frequency representation: 

G C pt(«w„) = G'{iuj n ) + G' '(iw„)VbGcpt(«w„) , (33) 

where ui n = (2n + 1)tt/ (3 with integer n, and fat symbols 
stand for matrices with respect to orbital indices only. 
After analytical continuation to arbitrary complex fre- 
quencies iuj n — > uj, we therewith recover the usual equi- 
librium CPT equatio n 22 ' 23 which may be solved by ma- 
trix inversion: 



G 



CPT 



G'H" 1 - V B 



(34) 



where translational symmetries of the lattice may be ex- 
ploited by Fourier transformation in addition. 

Eq. (|32|) holds for the equilibrium but also for the gen- 
eral non-equilibrium case, i.e. for time inhomogeneous 
Green's functions. Physically, it is a consequence of 
causality since the preparation of the initial state cannot 
depend on the subsequent time evolution of the system. 

The CPT does respect this condition: Consider an ex- 
pression of the form 

I(t,t')= f alt" f dt"'G 1 {t,t")'£{t",t"')G2(t"',t') , 



(35) 

as it occurs in the Dyson equation (|23[) or, in a simpler 
form, in the CPT equation (|3U|) . and assume the external 
time variables t and t' to lie on the Matsubara branch. 
After integrating over t'" , the integrand for the remaining 
t" integration depends on t" and t, t' only. In particular, 
since t, t' by assumption are always "later" than t" on 
7, if t" is real, it does not matter whether t" lies on the 
upper or on the lower branch of 7. Therefore, the integra- 
tion along the entire Keldysh branch does not contribute 
to the integral and 

I(t,t') = ^ dt" f dt"'G 1 (t,t")Il(t",t'")G2(t'",t') . 

(36) 
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Using the same arguments, we can then also replace 

\ dt!" ^ / dt" 1 (37) 

J 7 J to 

and we are left with integrations along the Matsubara 
branch only. 

B. Time discretization 

The numerical evaluation of the non-equilibrium CPT 
proceeds in two steps: (i) The contour-ordered Green's 
function G' of the reference system of disconnected clus- 
ters has to be calculated. If the individual cluster is 
sufficiently small, this can be done by full diagonaliza- 
tion of B and H (t). The computation is straightforward, 
(ii) The CPT equation ([30) must be solved. This is a 
Fredholm integral equation of the second kind which has 
the same formal structure as Dyson's equation ([2"3") . The 
standard approach consists in a discretization of the time 
variables to cast the CPT equation into a matrix form 
and to employ standard techniques for the solution of in- 
homogeneous linear systems of equations for its solution. 
It is recommcndablc to consider the CPT equation ([3TJ) 
in the form 

(1 - G' ■ V K ® 1)G CPT = G' , (38) 

as its solution formally requires a single inversion of a 
well-conditioned matrix only. 

We use ./Vm time slices for the Matsubara branch and 
jVk time slices for the upper as well as for the lower 
branch. This leads to a matrix dimension of Nm + 2Nk- 
Using Eq. ([5"2") to separate the solution of the CPT equa- 
tion on the Matsubara branch from the rest of the prob- 
lem, leads to Nm x Nm and 2Nk x 2Nk matrices only. Ex- 
ploiting further properties of the contour-ordered Green's 
function, one can reformulate Dyson's equation such that 
only Nm x Nm, Nk x Nk and Nm X TVk matrices must 
be considered for five independent quantities . 17 > 35 

For the time discretization, a maximal real time i max 
has to be introduced as a cutoff of the Kcldysh contour. 
This can be justified with arguments analogous to those 
given in the preceding section: If t, t' < i max , the inte- 
grations over t" (and t" 1 ) in Eq. (|3"5j) and thus in Eq. 
(|23p from t" = t max to t" = oo (upper branch) and from 
t" = oo to t" = t max (lower branch) cancel each other. 
Hence, any choice of i max > is justified. On the other 
hand, i max determines the maximal observation time up 
to which the Green's function GcpT,oa' (t, t') and thus ex- 
pectation values {A) t can be calculated. An immediate 
consequence of this is that non-equilibrium CPT can- 
not access the long-time behavior of observables: The 
numerical effort is dominated by the solution of linear 
systems of equations with a dimension Nk proportional 
to t ma x and therefore increases asymptotically as £j^ ax . 
Note, that matrix dimensions also increase due to site 
and orbital indices. 



C. Limiting cases 

Comparing the exact with the CPT self-energy, Eq. 
([25) with Eq. ([2D, shows that the non-equilibrium CPT 
becomes exact in the non-interacting limit Uk = as 
well as in the limit of decoupled clusters Vk = 0. The 
latter is, of course, trivial. The non-interacting limit, on 
the other hand, provides a serious check for the numerical 
evaluation of the theory. 

We have performed calculations for the single-band 
Hubbard model on a linear chain consisting of L sites 
with open boundaries: 

L-l u 
i=l c=t,4- 

(39) 

Here, T = 1 is the nearest-neighbor hopping which fixes 
the energy scale. Using the non-equilibrium CPT for 
Hubbard interaction U = 0, we have calculated the site- 
dependent occupation (n^t = (cj cr Cj cr )t = rii(t) as a 
function of the time t for spin-symmetric conditions. The 
initial state, prepared at t = to = 0, is assumed to be 
a pure state where N — L electrons occupy the sites 
i = 1, ...,L/2. This is a half-filled chain with all electrons 
located on the left half. Calculations are performed for 
L = 4 sites to allow for a check of the CPT results against 
the exact time evolution of rii(t). The reference system is 
taken to be given by two (non-interacting) Hubbard clus- 
ters consisting of two neighboring sites each such that the 
inter-cluster hopping, which in the CPT is treated pcr- 
turbativcly to all orders, is given by the hopping between 
the right site of the first and the left site of the second 
cluster. 

Fig. @] shows the results for different At = t mELX /NK- 
Choosing i max = 10 and At = 0.01 implies Nk = 2000 
time points on the Keldysh branch. As can be seen from 
the figure by comparing with the exact solution, this 
turns out to be sufficient for convergence of the results. 
The figure also demonstrates that the numerical evalua- 
tion recovers the {7 = limit correctly. The physics of 
this example is simple: For small t, the occupation of the 
second site quickly decreases, while due to Pauli block- 
ing, the occupancy at the first site starts to decrease with 
some time delay On a larger time scale, a strongly os- 
cillatory time evolution is observed as it is characteristic 
for a finite small system. 



D. Initial-state correlations 

For the above calculations, we only took the Keldysh 
contour into account and set f3 — 0. This is correct for an 
initial state represented by a Hamiltonian B with vanish- 
ing inter-cluster hopping Vb as it is the case here: The 
initial state is obtained as the ground state of the Hub- 
bard model with vanishing hopping between sites 2 and 
3 and suitably chosen on-site energies to realize a filled 
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FIG. 4: Time dependence of the site occupations in a four- 
site Hubbard chain for U — 0. Results of the non-equilibrium 
cluster-perturbation theory (lines) are compared with the ex- 
act result (points) for different time discretizations At on the 
Keldysh contour. The initial state is displayed schematically. 
At t = four electrons occupy the sites 1 and 2 while sites 
3 and 4 are empty. The CPT treats the nearest-neighbor 
hopping between sites 2 and 3 perturbatively to all orders. 



left and an empty right cluster. 

Vb = implies that the Matsubara branch is irrelevant 
for the time evolution within the CPT. To prove this, we 
consider the CPT equation Eq. (f3"U)) . As a matrix in t,t' 
the Green's function of the reference system consists of 
four blocks, 



G' = 



^KK 




"MK 





(40) 



where K refers to the upper and the lower branches of the 
Keldysh contour and M to the Matsubara branch. The 
block structure for non-retarded, instantaneous potential 
scattering is simple: 



V K 



V® 1 


0' 





o t 



(41) 



The matrix is diagonal and the MM block is zero for an 
initial state with Vb = 0. This immediately implies that 
the KK block of the CPT Green's function satisfies a 
simplified CPT equation, 

Gcpt,kk = G' KK + G' KK ■ V ® 1 • GcPT.KK ; (42) 

and depends on the KK block of the reference system's 
Green's function only. 

Within general non-equilibrium perturbation theory, 
the Matsubara branch cannot be disregarded unless the 
initial state is uncorrelated; 17 ' 35 Only if B\ = there are 
no vertices with imaginary time entries in the diagram- 
matic expansion of G. In the presence of initial-state 
correlations, however, the Matsubara branch is needed 
to expand the many-body density operator in terms of 
the non-interacting density operator which is a necessary 
prerequisite for the application of Wick's theorem. 



0.60 



u 0.55 - 



■& 0.50 



0.45 - 



0.40 




2 3 4 

time t [1 / n.n. hopping] 



FIG. 5: Time dependence of the spin-t site occupations in 
a four-site Hubbard chain for U — as obtained by non- 
equilibrium cluster-perturbation theory starting from a ref- 
erence system with decoupled two-site clusters and treating 
the nearest-neighbor hopping between sites 2 and 3 pertur- 
batively to all orders. The initial state is the ground state of 
the U — chain at half-filling with a local magnetic field of 
strength —h at site 1 and +h at site 2 with h — 0.2. CPT 
results for large j3 = 10 and for j3 = (i.e. neglecting the 
effects of the Matsubara branch). Exact results (points) are 
shown for comparison. 



Within non-equilibrium CPT, on the other hand, in- 
teraction vertices generated by Uk (including Ub) are 
taken into account to all orders for Vk = by the nu- 
merically exact calculation of the Green's function of the 
reference system G'. The subsequent summation of di- 
agrams generated by Vk, however, can be restricted to 
vertices on the Keldysh contour only since Vb = is as- 
sumed. The absence of effects of initial-state correlations 
on the real time evolution must therefore be seen as an 
artifact of the CPT. In fact, the self-energy diagram (c) 
in Fig. |3] is just a prime example to sec this: Wc assume 
the interaction vertices in this diagram to have imagi- 
nary time entries, i.e. we assume the interaction lines to 
be labeled by Ub, which may occur in case of a corre- 
lated initial state. Now, while the diagram is neglected 
within CPT, it gives a non- vanishing contribution within 
full perturbation theory even if Vb = since an interac- 
tion vertex at imaginary time and a potential-scattering 
vertex at real time can be connected by a non-vanishing 
clement of the MK block of the free propagator G' . 

Fig. [5] gives an example for a case where, within CPT, 
the effect of the Matsubara branch is essential. We again 
consider the U = Hubbard chain with L = 4 sites at 
half-filling. The system is assumed to be initially in the 
ground state of the same model but with an external 
magnetic field. This might also be seen as a magnetic- 
field quench. To induce a spatially asymmetric situation, 
we consider an additional field term 



H 



field 



= fe£(-i)V t 



(43) 
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to the Hamiltonian Eq. (|39|) which is staggered and non- 
zero on sites 1 and 2 only. In the figure, the resulting 
exact time dependence for t > to = is shown for {ni^) t 
as points. For the a =\. channel we have (n.jj,)t = 1 — 
(riif)t' For t = 0, the magnetic moments at sites 1 and 2 
are considerably larger than those at sites 3 and 4, due 
to the locally applied field. For t = 3, the situation is 
reversed, and the moments on sites 3 and 4 are larger. 

Since U = the CPT is expected to provide the exact 
result. In fact, for a reference system with decoupled 
two-site clusters (1 and 2 decoupled from 3 and 4) the 
calculation for j3 = 10 is close to the exact solution. This 
holds for the initial state, as can be seen be comparing 
the site occupations with the exact ones for t = 0, as well 
as for the subsequent time evolution. Residual deviations 
result from the finite time grid with At = 0.005 on the 
Kcldysh and on the Matsubara branch. 

In addition, the result of a CPT calculation with /3 = 
is shown in Fig. [5] This corresponds to a calculation 
on the Keldysh branch only but starting with the same 
Green's function of the reference system. Obviously, 
there are strong deviations from the exact result which 
proves the relevance of the Matsubara branch for the 
CPT calculations. For the present example the initial 
state is given as the ground state of a Hamiltonian B 
with Vb =V^0. Consequently, the Matsubara branch 
is required to restore the effect of the inter-cluster poten- 
tial in the initial state. Note that a finite field strength is 
necessary here (h = 0.2). For h — > oo the two clusters of 
the initial-state Hamiltonian B would decouple dynami- 
cally, and the initial state could be described with Vb = 
cquivalently, and the Matsubara branch would become 
irrelevant. Furthermore, we note that the /3 = results 
correspond to a calculation with Vb = in the initial- 
state Hamiltonian B since, as argued above, in that case 
the simplified CPT equation (|42|) on the Keldysh contour 
holds. 



to a mixed thermal initial state with inverse temperature 
(3 < oo. Hence, the limits (3 — > oo and t — s- —i/3 do not 
commute. Even for large but finite j3 in the CPT calcu- 
lation, however, the behavior of G' aa ,(t,t') for t —> —i/3 
or t' — > —i/3 cannot be neglected, provided that the Mat- 
subara branch is necessary at all, of course. The reason 
is that G' aa , (t, t'), considered as a matrix in t and t', does 
not adopt a block-diagonal structure in the j3 — > oo limit. 
We therefore enforce the symmetries Eq. (|44|) and Eq. 
([45]) by hand: The expectation value with \^f) is calcu- 
lated for t = —it and t' = —ir' with < r, r' < (3/2 and 
the symmetry relations are then used to get G aa > (t, t') for 
(3/2 < t,t' < (3. Clearly, for finite (3 this introduces arti- 
ficial discontinuities of the Green's function at t = — i/3/ 2 
and t' = —i/3/2. The height of the jumps, however, dis- 
appear asymptotically for f3 — > oo. Consequently, it is 
easily verified numerically that convergence to the exact 
result can be achieved for large (3 if the symmetries Eq. 
(|1H and Eq. (|45|) are enforced while strong deviations 
from the exact result remain present even for (3 — > oo 
otherwise. 

For efficiency reasons, one may exploit more symmetry 
relations. In fact, we find it convenient to profit from the 
exact relations 

G aa '(t v ,t' A ) = -G* a , a (t' v ,t A ) , (46) 
G aa i(t A ,t v ) = —Ga> a (t A ,t v ) , (47) 

which hold for t, t' on the Keldysh branch. Here, t A in- 
dicates that t belongs to the upper branch while t v = t 
but lies on the lower branch. We also make use of time 
homogeneity on the Matsubara branch, 

G aa ,(t,t')=G aa ,(t-t'), (48) 

valid for t = —it, t' = —ir' and < r, r' < (3. 



E. Exploiting symmetries 

If the time evolution of a pure state p — \^)(^>\ is con- 
sidered, the symmetries of the contour-ordered Green's 
function must be taken into account carefully. In the 
CPT calculation, the pure initial state \^} is obtained as 
the ground state of a suitably chosen Hamiltonian B by 
exact diagonalization. |W) is then used to get the refer- 
ence system (cluster) Green's function G' as an expecta- 
tion value. On the other hand, the Matsubara branch has 
to be cut off at a finite parameter (3 < oo. This implies 
that for t — > — i(3 the cluster Green's function, obtained 
as a ground-state expectation value, cannot respect the 
symmetry relations 

G aa ,(t,t') = G* a , a (t',-t(3-t) , (44) 
G aa ,(t',t) = G* a , a (-iP-t,t') , (45) 

which hold exactly for t on the Matsubara and t! on the 
Keldysh branch and for a Green's function corresponding 



F. Short-time dynamics 

Fig. |6] shows the results of a calculation for the Hub- 
bard model Eq. (|39[) in the strong-coupling regime for 
U = 8. To allow for a comparison of the results from 
non-equilibrium CPT with the exact results, we consider 
the L = 4 site chain again. Initially, the system is pre- 
pared in the Neel state \^) where (^\ni^\^) = 1 for i = 1 
and i = 3 and (^\rii^\^i) = for i = 2 and i = 4 and 
where (*|n»4,|*) = 1 - {^\m^). For strong U at half- 
filling the Hubbard model maps onto the antiferromag- 
netic Heisenbcrg model with a ground state and excited 
energy eigenstates that are different from the classical 
Neel state. This induces a non-trivial dynamics as can 
be seen from the exact calculation (blue lines) in Fig. [5] 

The Neel state may be obtained as the ground state of 
an initial-state Hamiltonian B with a staggered magnetic 
field term as in Eq. (|4"3"f but applied to all sites and with 
field strength h — > oo. This implies that the sites are 
decoupled dynamically, and that Vb = can be assumed 
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FIG. 6: Time evolution as obtained from non-equilibrium 
CPT (red), exact time evolution (blue) and time evolution 
of the reference system (green) in the strong-coupling regime 
(U — 8) of a half-filled four-site Hubbard chain with a Neel 
initial state. Reference system: decoupled two-site clusters. 
Inter-cluster potential: nearest-neighbor hopping connecting 
sites 2 and 3. 



G. Coupling to an infinite bath 

The (equilibrium) CPT has actually been designed to 
treat correlated electrons on an infinite lattice. For the 
non-equilibrium case, the results presented above repre- 
sent simple test calculations which demonstrate that two 
correlated clusters with L c sites each can be coupled to 
a single but larger cluster with 2L C sites. This scheme 
can be iterated straightforwardly to build up extended 
lattices with or without translational symmetries. 

Besides this, the non-equilibrium CPT can also be used 
to couple a small correlated "system" to an uncorrelated 
"bath" with a large number of degrees of freedom, such 
as a magnetic nanostructure on a metal surface or a 
molecule coupled to external leads etc. We assume that 
the Hilbcrt-space dimension of the correlated system is 
sufficiently small such that the contour-ordered Green's 
function G' s can be calculated exactly. By means of Eq. 
([24| and Eq. (|25j) we also have the Green's function of 
the bath G' B for an in principle arbitrarily large number 
of uncorrelated bath sites. Hence, the Green's function 
of the decoupled reference system, given by B' and H'(t), 
can be written as a 2 x 2 matrix 



for the initial state. Consequently, we are allowed to 
disregard the Matsubara branch. 

For the CPT calculation we again start from a refer- 
ence system with decoupled two-site clusters. By con- 
struction, the initial state is described correctly within 
the CPT approach. As can be seen from Fig. |6l the site 
occupations obtained by CPT (red lines) deviate from 
the exact results for t > as expected for U ^ 0. For 
comparison, the time dependence of the site occupations 
of the reference system are given in addition (green line) . 
The reference system has a higher symmetry which leads 
to occupations of sites 1 and 2 that are related to each 
other by spin reversal. For larger times the CPT results 
seem to follow more or less the time dependence of the 
site occupations in the reference system. This means that 
the approximation is not able to describe the effects of 
inter-cluster correlations correctly and that intra-cluster 
effects dominate the behavior of {ni a )t at large times. 

On the other hand, at short times t < 0.3, the CPT 
results arc clearly different from the site occupations of 
the reference system and to a high precision follow the 
exact trend. We conclude that inter-cluster correlations, 
as represented by the diagram (c) in Fig. [3l are ineffec- 
tive at short times even if the interaction is strong. The 
fact that the non-equilibrium CPT describes the short- 
time dynamics of single-particle operators exactly, is in- 
terpreted to be the analog of the fact that the equilibrium 
CPT predicts global, i.e. frequency-integrated properties 
of the single-particle excitation spectrum correctly. The 
CPT apparently respects to a good approximation the 
first few non-equilibrium moment sum rules which deter- 
mine the short-time dynamics ! 36 ! 37 



G' 











G'b 



(49) 



with entries referring to system or bath orbitals. The 
coupling of the system to the bath is provided by the 
"inter-cluster" term, i.e. by the hybridization 



V K 






V K ' 








(50) 



Within non-equilibrium CPT, the Green's function of the 
entire system, 



G 



CPT 



Gcpt,s 


CcPT,SB 


GcPT.BS 


GcPT,B 



(51) 



is obtained from the general CPT equation, Gcpt = 
G' + G' ■ V K (8 1 • Gcpt (see Eq. ) . With the definition 
of the non-equilibrium hybridization function, A = Vk §5 
1 • G' B (Vk (8 l) f , or, in short, 



v ■ Gr ■ 



(52) 



the CPT Green's function of the system's degrees of free- 
dom is obtained as: 



Gcpt,: 



(53) 



For the situation considered here, this can be seen as 
a simplified CPT equation which is decoupled from the 
remaining CPT equations for the bath and system/bath 
Green's functions, 



V ■ Gcpt,: 



G' c 



(54) 
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and 

V ■ GcPT,BS = A • GcPT,S ) 

Gcpt.sb • V* = G' s ■ (V • G C pt,b • Vt) . (55) 

Eqs. (|53|) - (|55[) have simple diagrammatic representa- 
tions known from scattering theory. 

For the following a small system with L c sites is con- 
sidered and a hybridization that links the site iq of the 
system to a single site io,bath of the bath. Let po(e) be 
the non-interacting local density of states of the bath 
at io.bath- This implies that the hybridization function 
A is non-zero at io only, and A(t,t') = A» j (i,t') = 
V(t)G' Bloc (t,t')V(t') where the local bath Green's func- 
tion at i ,bath is given by 

f e -i{e-n)(t-t') 

iG' B , loc .(t,t') = J fePo{e) 1 + e _ fKe _ ll) (56) 
if t is "later" then t' on the contour, and 

f e -»(e-|*)(*-t') 

iG' B , loc Xt,t ') = -J depo(e) eP(e _ /l) + l (57) 

if t' is "later" then t on the contour. This means that the 
bath is fully characterized by its local density of states 
Pq(s) at «o,bath- The CPT equation Eq. ([53]) then pro- 
vides the system's Green's function at 

GcPT,S,i io = 7^=1 T" ) (58) 

where fat quantities are matrices in t, t' only. For the 
other sites we have: 

G C PT,s,ij = G'ax + G' Stiio A— — -^G' s - loi . (59) 

S.ioio 

For the numerical calculations we consider a system 
in a linear geometry with L c > 2 sites. The Hubbard 
interaction U is non-zero at sites i = 1 and i = 2 only, and 
the hopping between nearest neighbors is set to T = 1 
to fix the energy and time scales. System sizes range 
from L c = 2 to L c = 8. The latter is the maximum 
size that can conveniently be treated by means of exact 
diagonalization. Via non-equilibrium CPT this system is 
coupled at the site i = io = L c to a bath with a semi- 
elliptic density of states po(s) of bandwidth W = AT. 
This is just the local density of states at the first site 
*o,bath for a semi-infinite linear chain. Both, the system 
and the bath, are taken to be at half-filling, i.e. we set the 
chemical potential p = 0, and assume vanishing on-site 
energies for all sites except for i = 1, 2 where the on-site 
energy is — U/2. In the ground state n„ = (m a ) =0.5 
for system and bath. 

However, the initial state is taken to be the ground 
state of another Hamiltonian B which differs from H by 
(i) the hopping between sites i = 2 and i — 3. This hop- 
ping is suddenly switched on at time t = 0. Furthermore, 
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FIG. 7: Pictorial representation of the final-state Hamilto- 
nian H governing the time evolution and the Hamiltonian B 
generating the initial state as its ground state. The system 
size is io — L c . The bath is a semi-infinite linear chain starting 
at io,bath = L c + 1. The coupling between system and bath V 
is treated within CPT (V = T). Blue circles represent sites 
with U > while open red circles stand for sites with U = 0. 
Initially (t — 0) there is a ground state of a an isolated two- 
site Hubbard cluster perturbed by a staggered field h coupling 
to spin or charge degrees of freedom. The field is suddenly 
switched off, and the coupling of the two-site cluster to the 
rest of the system is switched on. Due to particle-hole sym- 
metry, the entire system ( "system" plus bath) is and stays at 
half-filling. 

(ii) the correlated two-site model for the initial state is 
perturbed by either a spin or by a charge excitation. This 
is realized by applying a respective staggered field term: 

B = B(0) - h sphl (n n - nil) + ^s P in("2t - n 2i) (60) 

or 

B = S(0)-/l cha rgc("lt+"-u)+ /l chargc("2t+"-24.) • (61) 

The Hamiltonians of the initial ground state and of the 
transient final state, i.e. B and H, are shown schemati- 
cally in Fig. [3 Note that the CPT describes the initial 
state exactly because the correlated sites are decoupled 
and because the coupling of the rest of the sites of the 
system to the bath is taken into account exactly via CPT 
since these sites are uncorrelated. Converged results are 
obtained with the choice /? = 5 for the Matsubara branch. 

Physically, we expect that the initial local perturbation 
propagates through the chain and dissipates into the bath 
such that the system relaxes to its ground state with 
n a = 0.5. For U = 2 Fig. [5] displays the result of a 
calculation for a spin excitation with /i S pin = 0.1, and 
Fig. [9] the results for a charge excitation with /i c harge = 
0.1. We find that the results for the time dependence 
of rif(t) = 1 — ni(t) improve with increasing size of the 
system. Clearly, if the CPT was exact there should not 
be any differences between the results of calculations for 
different L c . 

At site i = 1 (Fig. [51 upper panel) the result obtained 
for the smallest system with L c = 2 shows a strongly os- 
cillating trend with hardly any damping despite the pres- 
ence of the bath. Here, the non-equilibrium CPT appears 
reliable on a short time scale t < 2 only as can be seen by 
comparing with L c > 4. By comparing with the result for 
the largest size, it is obvious that this time scale rapidly 
grows with increasing system size. For L c = 8 and up 
to the accessible maximal time i max = 10, the trend of 
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FIG. 8: Time evolution of the average occupation rif — 
{nit) = 1 — ( n il) a t site i — 1 (upper) and site i — L c 
(lower panel) of the system displayed in Fig. after a spin 
excitation using a staggered field h S pm = 0.1 in the initial- 
state Hamiltonian (see Eq. (|60|) V Insets: Comparison of the 
CPT result for L c = 8 with the corresponding result obtained 
from a calculation for an isolated system without bath. The 
CPT calculations have been performed for different system 
size L c as indicated. Further parameters: [7 = 2, half-filling, 
T — V — 1. 



n-f(t) follows our expectation: The initial spin polariza- 
tion quickly decreases and, apart from weak remaining 
oscillations, approaches — n\. ~ 0. 

The lower panel of Fig. [5] shows n-j-(i) at site i = in = 
L c . As the local spin excitation requires a finite time to 
propagate to i — iq, the results for the different system 
sizes show a response that is more and more delayed with 
increasing L c . The excitation is nicely seen to propagate 
through in and leaving the site in an almost unperturbed 
state thereafter. The upturn of n-j-(t) for t > 7 in the 
calculation for L c = 4, however, must be interpreted as 
an artifact. Here the system size is still insufficient to 
predict the correct trend up to t ma x- On the other hand, 
one should note that the decrease in the amplitude of 
the response with increasing L Cl and thus with increasing 
distance from the initial perturbation, is reasonable. 

Fig. [HI presents the time-evolution of n(t) = n^{t) = 
n^(t) after an initial charge excitation at i = 1. Again, 
the initially strong deviation from the equilibrium value 
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Time evolution of the average occupation n = 
(riii) at site i = 1 of the system displayed in Fig. 
[7] after a charge excitation with /icharge = 0.1 (see Eq. (|60|l 'l. 
System size L c as indicated, U = 2, half-filling, T = V = 1. 



n(t) = 0.5 is quickly dissipated to the bath while the re- 
maining low-amplitude oscillations are expected to decay 
on a time scale beyond i max . We also find that the CPT 
results rapidly improve with increasing system size L c . 

Comparing the results for spin and charge excitations, 
we note that the system is substantially more susceptible 
to a staggered field that couples to the spin as compared 
to a field coupling to the charge degrees of freedom; for 
^-charge = ^-spin we find oscillations with larger amplitudes 
in Fig. |8j Furthermore, the characteristic frequency of 
the oscillations seen in Fig. [3] for the spin excitation is 
clearly smaller than the corresponding one for the charge 
excitation (Fig. HJ). These facts are strongly dependent 
on U. With increasing U we find that the characteris- 
tic frequency for the spin excitation is roughly given by 
w S pin ~ 1/f which corresponds to the low-energy Hcisen- 
berg scale, while for the charge excitation w c haxge ~ U 
which corresponds to the high-energy Hubbard bands. 
This is accompanied by an increase (decrease) of the am- 
plitudes for the oscillations following the spin (charge) 
excitation. For strong U, the system is very weakly sus- 
ceptible to a perturbation coupling to the charge as com- 
pared to the spin degrees of freedom. 

We conclude that the non-equilibrium CPT is in fact 
able to describe the dynamics following a perturbation of 
a small correlated system in a non-interacting environ- 
ment and the dissipation of a local spin or charge excita- 
tion to a large uncorrelatcd bath. It is important to note, 
however, that the above-mentioned effects are to some 
extent already captured in a calculation for V — in the 
final state, i.e. in a calculation without bath. This is most 
apparently seen in the inset of the upper panel in Fig. [SI 
where rif(t) obtained by CPT is compared with the re- 
sult for the isolated cluster at L c = 8. The CPT does 
improve the calculation for the isolated cluster but the 
gain is small. The reason is that the "reflection" of the 
propagating excitation at the boundary i = ig and the 
back-propagation to i = 1 takes a time close to t — i roax . 
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On the other hand, at site i — iq (see lower panel), the 
CPT substantially improves the isolated-cluster calcula- 
tion by predicting a much stronger damping. 

These observations can also be understood from the 
diagrammatic construction of the non-equilibrium CPT 
(see Fig. |3]) by assuming that non-diagonal elements of 
the free intra-cluster propagator, G ^ with i =/= j but 
i,j G / are small compared to diagonal elements G' u 
and decrease with increasing distance \i— j | . The diagram 
to the self-energy in Fig. neglected within the CPT, 
necessarily involves two non-diagonal propagators G' 4 • 
with i = 1 or i = 2 and j = iq since the U and the V 
vertex are local and separated by a distance io — 1 (see 
Fig- Ell- It is therefore of the order O(G' iig ) and vanishes 
with io — > oo. The same argument can be given for any 
vertex-correction diagram and hence the CPT becomes 
exact in the limit of L c — > oo, as expected. 

Likewise, we can argue that the contribution of ne- 
glected vertex-correction diagrams to the site occupation 
at i = 1 or i = 2 are of the order O(G' 0iia ). On the 
other hand, for i = io, the CPT provides a more reli- 
able estimate since vertex corrections are already of the 
order 0(G'q u ) because of the necessary two additional 
non-diagonal propagators. 

V. CONCLUSIONS AND OUTLOOK 

Usually, the cluster-perturbation theory is seen as the 
first non-trivial level in a systematic strong-coupling ex- 
pansion, i.e. an expansion in the inter-cluster hopping V 
around a state with decoupled clusters but finite and ar- 
bitrarily strong Hubbard-type interaction U. Here, we 
have shown that the same CPT can be recovered strictly 
within the framework of weak-coupling perturbation the- 
ory. This is achieved by formally summing up certain 
classes of diagrams that are generated when treating V 
and U pcrturbatively. In this way the CPT can be inter- 
preted as an approximation that neglects vertex correc- 
tions, i.e. the influence of scattering at the inter-cluster 
potential on the renormalization of propagators due to 
the interaction. One of the benefits of this reformulation 
is that therewith one can straightforwardly extend the 
CPT to study the real-time dynamics of systems far from 
equilibrium. One simply has to copy the formalism and 
paste it to the Keldysh-Matsubara time contour. This 
defines the non-equilibrium CPT studied here. 

The non-equilibrium CPT is characterized as follows: 
(i) It comprises the conventional CPT for the description 
of the initial thermal state and fully reduces to conven- 
tional CPT in the case of thermal equilibrium, i.e. for 
the case where the Hamiltonian H(t) that determines 
the time evolution is assumed to be time independent 
and set equal to the Hamiltonian B that defines the ini- 
tial thermal state. 

(ii) The non-equilibrium CPT respects the physical 
consequences of the causality principle: Within the CPT 
the time evolution of the system depends on the initial 



state preparation but not vice versa. 

(iii) The approach is rather flexible and can be applied 
to a large class of models, namely lattice fermion (or 
boson) models with local Hubbard-type interactions in- 
cluding impurity models such as the single-impurity An- 
derson model. For bosons, however, the treatment of the 
condensate phase requires additional efforts.— Further- 
more, systems with non-local interactions, like a nearest- 
neighbor density interaction cannot be treated without 
further approximations, such as a mean-field decoupling 
of inter-cluster interaction terms.— 

(iv) Due to the necessity to solve a generalized CPT 
equation for time-inhomogencous Green's functions, op- 
erations involving objects indexed with two discrctized 
time variables have to be performed. This limits the nu- 
merical evaluation of the scheme to short and interme- 
diate time scales in practice. On the other hand, there 
are in principle no limitations concerning the time depen- 
dence of the Hamiltonian, and the non-equilibrium CPT 
can likewise treat sudden parameter changes or periodi- 
cally driven systems, for example. 

(v) The neglect of vertex corrections represents a se- 
vere approximation. This approximation is in principle 
controlled, however, by the cluster size, i.e. the (non- 
equilibrium) CPT approximation improves with increas- 
ing L c . This is shared with the conventional (thermal) 
CPT and classifies the scheme as a cluster mean-field ap- 
proach where correlations are treated exactly up to the 
cluster extension and treated in a mean-field way be- 
yond this scale. For impurity-type models with a single 
or a few correlated sites and a continuum of uncorrected 
bath degrees of freedom, the approximation has also been 
seen to improve with increasing distance of the correlated 
sites from the cluster boundary. Here, "improvement" 
means that the dynamics of expectation values of single- 
particle observables can be traced reliably on longer and 
longer time scales. On a very short time scale, the non- 
equilibrium CPT has been found to recover the exact 
solution, i.e. it apparently (like the equilibrium CPT) re- 
spects the first non-trivial moment sum rules. 

(vi) The non-equilibrium CPT can also be character- 
ized as a scheme that interpolates between the isolated- 
cluster limit (V = 0) and the band limit (U = 0) which 
are recovered exactly. However, already at the second or- 
der in the interaction strength there are diagrams miss- 
ing. An interesting case that should be accessible to the 
method and has been studied experimentally, for bosonic 
atoms in optical lattices ; 40 ' 41 are weakly coupled double 
wells or weakly coupled plaquettcs. 

Concluding, the approach represents a very flexible 
and easy to handle method with very moderate com- 
putational cost that can give a first access to a rather 
broad class of systems of strongly correlated electrons 
far from equilibrium. On the other hand, its main draw- 
backs consist in the missing self-consistency, the neglect 
of correlations beyond the cluster size and also the arti- 
ficial breaking of lattice symmetries. The present work 
has presented a number of test calculations. These can 
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be improved in various ways to achieve more reliable re- 
sults: Larger clusters can be taken into account by re- 
placing the exact-diagonalization approach for the com- 
putation of the Keldysh Green's function with e.g. a time- 
adaptive Krylov construction^ A (strong-coupling) dia- 
grammatic expansion around the non-equilibrium CPT 
may be used^i to include some of the neglected vertex 
corrections. Alternatively, one can also attempt to en- 
large the class of diagrams considered in the presented 
weak-coupling expansion. Finally, an optimization of 
intra-cluster one-particle parameters can be envisaged to 
introduce a self-consistent feedback within the method 
which is necessary to study phase transitions and to 
make contact with non-equilibrium dynamical mean-field 



theory, for example. This can be accomplished by a 
suitable generalization of the sclf-energy-functional ap- 
proach. Work along these lines is in progress.-^ 
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